Done.
[ NOTE, FIGURE ]: Red horizontal line indicates a signed R^2 of 0.9
Setting soft-thresholding power to: 15.
Power-transforming the gene-gene similarity matrix...Done.
---------------------------------------------------
4. Convert into topological overlap matrix (dissTOM)
---------------------------------------------------
Creating dissTOM...Done.
Performing hierarchical clustering on dissTOM...Done.
---------------------------------------------------
5. Identify modules (clusters)
---------------------------------------------------
Merging modules that have a correlation ≥ 0.9 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.
Module (cluster) size:
mergedColors
black cyan darkgreen darkgrey darkorange
1392 158 88 81 78
darkred darkturquoise greenyellow grey60 lightgreen
105 87 204 133 127
lightyellow magenta midnightblue orange paleturquoise
127 211 150 80 56
pink purple red royalblue saddlebrown
212 204 234 124 59
skyblue steelblue tan turquoise violet
72 57 198 320 54
white yellow
77 255
Merging modules that have a correlation ≥ 0.85 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.
Module (cluster) size:
mergedColors
black cyan darkgreen darkgrey darkred
1392 236 88 605 105
darkturquoise grey60 lightgreen magenta orange
214 283 438 211 80
paleturquoise pink royalblue saddlebrown skyblue
56 416 124 59 72
steelblue tan violet yellow
57 198 54 255
Merging modules that have a correlation ≥ 0.8 ...Done.
[ NOTE, FIGURE ] Plotting identified clusters before and after merging.
Module (cluster) size:
mergedColors
black darkgreen darkgrey darkred darkturquoise
1883 88 605 105 214
lightgreen magenta orange royalblue saddlebrown
1108 211 80 407 59
skyblue steelblue violet
72 57 54
Cutoff used: 0.8
Number of modules identified: 13
Calculating module-module similarity based
on module-eigengene-expression...Done.
Tidying module names...Done.
Plotting adjacency matrix for module-module similarity...
Visualizing a simplified representation of the circadian GCN
Annotate the network
Identify rhythmic modules
Show code
# Obtain a list of genes in each modulel_module_genes <- mods[["module_genes"]] |>arrange(module_identity) |>group_split(module_identity) |> purrr::map(~ .x |>pull(gene_name) ) |>setNames(unique(mods[["module_genes"]]$module_identity))# Load results of rhythmicity analysesdb_rhy <- purrr::map( sample.names,~load_rhy_genes(sample = .x )) |>setNames(sample.names)# Set your p-value of choice###-###-###-###-###-###-###-###-# col_pval = "BH.Q"col_pval ="default.pvalue"# col_pval = "raw.pvalue"###-###-###-###-###-###-###-###-# Obtain a list of rhythmic genes in each tissuel_rhy_genes <- purrr::map( sample.names,~ db_rhy[[.x]] |> purrr::map(~ .x |>filter(if_all(all_of(col_pval),~ .x <0.05 ) ) |>filter( ID %in%unlist(l_module_genes) ) |>pull(1) |>unique() ) |> purrr::compact()) |>setNames(sample.names)
Modules vs. Rhythmic genes
Show code
# LIST ONE - WGCNA moduleslist1 <- l_module_genessapply(list1, length) |>print()
“We begin by calculating the intramodular connectivity for each gene. (In network literature, connectivity is often referred to as”degree”.) The function intramodularConnectivity computes:
the whole network connectivity kTotal,
the within (intra)module connectivity kWithin,
the extra-modular connectivity kOut=kTotal-kWithin, and
the difference of the intra- and extra-modular connectivities kDiff = kIn - kOut = 2*kIN-kTotal
Show code
# From what I can tell, colorh1 in the tutorial refers to moduleColorscolorh1 <- mods[["modules"]]$colorsadj_matrix <- mods[["adj_matrix"]]# Calculate the connectivities of the genesAlldegrees1 = WGCNA::intramodularConnectivity(adjMat = adj_matrix,colors = colorh1) |>mutate(gene_name =rownames(adj_matrix),across(matches("^k"),~round(.x, 2) ) ) |>glimpse()
# tidydata.db |> # purrr::map_int(# ~ nrow(.x)# )# We work with two sets:nSets =2;# Find the common set of genes across all samplescommon_genes <- tidydata.db |> purrr::map(~ .x |>pull(gene_name) |>unique() ) |> purrr::reduce( intersect ) |>intersect( mods[["module_genes"]]$gene_name ) # filter to keep only genes used in creating REF networksamples.list <- tidydata.db |> purrr::map(~ .x |>filter( gene_name %in% common_genes ) )test_samples <-setdiff(sample.names, which.sample)l_mp <- purrr::map( test_samples,function (x) { multiExpr <-prep_data_module_preservation(data = samples.list,ref = which.sample, test = x )# CHECKS#---#---#---#---#---#---#---#---#---#---#---#---# Check that the data has the correct format for many functions # operating on multiple sets: exprSize = WGCNA::checkSets(multiExpr)# Check that all genes and samples have sufficiently low numbers of # missing values. gsg = WGCNA::goodSamplesGenesMS(multiExpr, verbose =1);cat("All samples okay?\n", gsg$allOK)#---#---#---#---#---#---#---#---#---#---#---#---# prep data#---#---#---#---#---#---#---#---#---#---#---#---# Expression data multiExpr_1 =list(ref =list(data = multiExpr[[1]]$data), test =list(data = multiExpr[[2]]$data) )## filter to keep only the genes that we are working with mod.identity <- mods[["module_genes"]] |>select( gene_name, old_labels ) |>filter( gene_name %in% common_genes ) |># !!this step is necessary!! #arrange(gene_name)## specify the module identity of the genes moduleColors <- mod.identity |>pull(old_labels) multiColor =list(ref = moduleColors)# Run module preservation#---#---#---#---#---#---#---#---#---#---#---#--- mp = WGCNA::modulePreservation( multiExpr_1, multiColor,referenceNetworks =1,nPermutations =200,calculateQvalue =FALSE,quickCor =0,verbose =1 ) mp })
Flagging genes and samples with too many missing values...
..step 1
All samples okay?
TRUE ..checking data for excessive amounts of missing data..
..calculating observed preservation values
..calculating permutation Z scores
..Working with set 1 as reference set
Flagging genes and samples with too many missing values...
..step 1
All samples okay?
TRUE ..checking data for excessive amounts of missing data..
..calculating observed preservation values
..calculating permutation Z scores
..Working with set 1 as reference set
Flagging genes and samples with too many missing values...
..step 1
All samples okay?
TRUE ..checking data for excessive amounts of missing data..
..calculating observed preservation values
..calculating permutation Z scores
..Working with set 1 as reference set
Flagging genes and samples with too many missing values...
..step 1
All samples okay?
TRUE ..checking data for excessive amounts of missing data..
..calculating observed preservation values
..calculating permutation Z scores
..Working with set 1 as reference set